Cut file at lines of the form "-- filename.ext --" and write to file
with that name.

-- ACTAN.FOR --
\begin{verbatim}
      FUNCTION  ACTAN(SINX,COSX)
      COMMON/C2/DE2RA,PI,PIO2,TWOPI,X3PIO2
      ACTAN=0.
      IF (COSX.EQ.0.  ) GO TO 5
      IF (COSX.GT.0.  ) GO TO 1
      ACTAN=PI
      GO TO 7
    1 IF (SINX.EQ.0.  ) GO TO 8
      IF (SINX.GT.0.  ) GO TO 7
      ACTAN=TWOPI
      GO TO 7
    5 IF (SINX.EQ.0.  ) GO TO 8
       IF (SINX.GT.0.  ) GO TO 6
      ACTAN=X3PIO2
      GO TO 8
    6 ACTAN=PIO2
      GO TO 8
    7 TEMP=SINX/COSX
      ACTAN=ACTAN+ATAN(TEMP)
    8 RETURN
      END
\end{verbatim}
-- DEEP.FOR --
\begin{verbatim}
*      DEEP SPACE                                         31 OCT 80
      SUBROUTINE DEEP
      COMMON/E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,
     1           XNDD6O,BSTAR,X,Y,Z,XDOT,YDOT,ZDOT,EPOCH,DS50
      COMMON/C1/CK2,CK4,E6A,QOMS2T,S,TOTHRD,
     1           XJ3,XKE,XKMPER,XMNPDA,AE
      COMMON/C2/DE2RA,PI,PIO2,TWOPI,X3PIO2
      DOUBLE PRECISION EPOCH, DS50
      DOUBLE PRECISION
     *     DAY,PREEP,XNODCE,ATIME,DELT,SAVTSN,STEP2,STEPN,STEPP
      DATA              ZNS,           C1SS,          ZES/
     A                  1.19459E-5,    2.9864797E-6, .01675/
      DATA              ZNL,           C1L,           ZEL/
     A                  1.5835218E-4,  4.7968065E-7,  .05490/
      DATA              ZCOSIS,        ZSINIS,        ZSINGS/
     A                  .91744867,     .39785416,     -.98088458/
      DATA              ZCOSGS,        ZCOSHS,        ZSINHS/
     A                  .1945905,      1.0,           0.0/
      DATA Q22,Q31,Q33/1.7891679E-6,2.1460748E-6,2.2123015E-7/
      DATA G22,G32/5.7686396,0.95240898/
      DATA G44,G52/1.8014998,1.0508330/
      DATA G54/4.4108898/
      DATA ROOT22,ROOT32/1.7891679E-6,3.7393792E-7/
      DATA ROOT44,ROOT52/7.3636953E-9,1.1428639E-7/
      DATA ROOT54/2.1765803E-9/
      DATA THDT/4.3752691E-3/

*     ENTRANCE FOR DEEP SPACE INITIALIZATION

      ENTRY DPINIT(EQSQ,SINIQ,COSIQ,RTEQSQ,AO,COSQ2,SINOMO,COSOMO,
     1         BSQ,XLLDOT,OMGDT,XNODOT,XNODP)
      THGR=THETAG(EPOCH)
      EQ = EO
      XNQ = XNODP
      AQNV = 1./AO
      XQNCL = XINCL
      XMAO=XMO
      XPIDOT=OMGDT+XNODOT
      SINQ = SIN(XNODEO)
      COSQ = COS(XNODEO)
      OMEGAQ = OMEGAO

*     INITIALIZE LUNAR SOLAR TERMS

    5 DAY=DS50+18261.5D0
      IF (DAY.EQ.PREEP)    GO TO 10
      PREEP = DAY
      XNODCE=4.5236020-9.2422029E-4*DAY
      STEM=DSIN (XNODCE)
      CTEM=DCOS (XNODCE)
      ZCOSIL=.91375164-.03568096*CTEM
      ZSINIL=SQRT (1.-ZCOSIL*ZCOSIL)
      ZSINHL= .089683511*STEM/ZSINIL
      ZCOSHL=SQRT (1.-ZSINHL*ZSINHL)
      C=4.7199672+.22997150*DAY
      GAM=5.8351514+.0019443680*DAY
      ZMOL = FMOD2P(C-GAM)
      ZX= .39785416*STEM/ZSINIL
      ZY= ZCOSHL*CTEM+0.91744867*ZSINHL*STEM
      ZX=ACTAN(ZX,ZY)
      ZX=GAM+ZX-XNODCE
      ZCOSGL=COS (ZX)
      ZSINGL=SIN (ZX)
      ZMOS=6.2565837D0+.017201977D0*DAY
      ZMOS=FMOD2P(ZMOS)

*     DO SOLAR TERMS

   10 LS = 0
      SAVTSN=1.D20
      ZCOSG=ZCOSGS
      ZSING=ZSINGS
      ZCOSI=ZCOSIS
      ZSINI=ZSINIS
      ZCOSH=COSQ
      ZSINH=SINQ
      CC=C1SS
      ZN=ZNS
      ZE=ZES
      ZMO=ZMOS
      XNOI=1./XNQ
      ASSIGN 30 TO LS
   20 A1=ZCOSG*ZCOSH+ZSING*ZCOSI*ZSINH
      A3=-ZSING*ZCOSH+ZCOSG*ZCOSI*ZSINH
      A7=-ZCOSG*ZSINH+ZSING*ZCOSI*ZCOSH
      A8=ZSING*ZSINI
      A9=ZSING*ZSINH+ZCOSG*ZCOSI*ZCOSH
      A10=ZCOSG*ZSINI
      A2= COSIQ*A7+ SINIQ*A8
      A4= COSIQ*A9+ SINIQ*A10
      A5=- SINIQ*A7+ COSIQ*A8
      A6=- SINIQ*A9+ COSIQ*A10
C
      X1=A1*COSOMO+A2*SINOMO
      X2=A3*COSOMO+A4*SINOMO
      X3=-A1*SINOMO+A2*COSOMO
      X4=-A3*SINOMO+A4*COSOMO
      X5=A5*SINOMO
      X6=A6*SINOMO
      X7=A5*COSOMO
      X8=A6*COSOMO
C
      Z31=12.*X1*X1-3.*X3*X3
      Z32=24.*X1*X2-6.*X3*X4
      Z33=12.*X2*X2-3.*X4*X4
      Z1=3.*(A1*A1+A2*A2)+Z31*EQSQ
      Z2=6.*(A1*A3+A2*A4)+Z32*EQSQ
      Z3=3.*(A3*A3+A4*A4)+Z33*EQSQ
      Z11=-6.*A1*A5+EQSQ *(-24.*X1*X7-6.*X3*X5)
      Z12=-6.*(A1*A6+A3*A5)+EQSQ *(-24.*(X2*X7+X1*X8)-6.*(X3*X6+X4*X5))
      Z13=-6.*A3*A6+EQSQ *(-24.*X2*X8-6.*X4*X6)
      Z21=6.*A2*A5+EQSQ *(24.*X1*X5-6.*X3*X7)
      Z22=6.*(A4*A5+A2*A6)+EQSQ *(24.*(X2*X5+X1*X6)-6.*(X4*X7+X3*X8))
      Z23=6.*A4*A6+EQSQ *(24.*X2*X6-6.*X4*X8)
      Z1=Z1+Z1+BSQ*Z31
      Z2=Z2+Z2+BSQ*Z32
      Z3=Z3+Z3+BSQ*Z33
      S3=CC*XNOI
      S2=-.5*S3/RTEQSQ
      S4=S3*RTEQSQ
      S1=-15.*EQ*S4
      S5=X1*X3+X2*X4
      S6=X2*X3+X1*X4
      S7=X2*X4-X1*X3
      SE=S1*ZN*S5
      SI=S2*ZN*(Z11+Z13)
      SL=-ZN*S3*(Z1+Z3-14.-6.*EQSQ)
      SGH=S4*ZN*(Z31+Z33-6.)
      SH=-ZN*S2*(Z21+Z23)
      IF(XQNCL.LT.5.2359877E-2) SH=0.0
      EE2=2.*S1*S6
      E3=2.*S1*S7
      XI2=2.*S2*Z12
      XI3=2.*S2*(Z13-Z11)
      XL2=-2.*S3*Z2
      XL3=-2.*S3*(Z3-Z1)
      XL4=-2.*S3*(-21.-9.*EQSQ)*ZE
      XGH2=2.*S4*Z32
      XGH3=2.*S4*(Z33-Z31)
      XGH4=-18.*S4*ZE
      XH2=-2.*S2*Z22
      XH3=-2.*S2*(Z23-Z21)
      GO TO LS

*     DO LUNAR TERMS

   30 SSE = SE
      SSI=SI
      SSL=SL
      SSH=SH/SINIQ
      SSG=SGH-COSIQ*SSH
      SE2=EE2
      SI2=XI2
      SL2=XL2
      SGH2=XGH2
      SH2=XH2
      SE3=E3
      SI3=XI3
      SL3=XL3
      SGH3=XGH3
      SH3=XH3
      SL4=XL4
      SGH4=XGH4
      LS=1
      ZCOSG=ZCOSGL
      ZSING=ZSINGL
      ZCOSI=ZCOSIL
      ZSINI=ZSINIL
      ZCOSH=ZCOSHL*COSQ+ZSINHL*SINQ
      ZSINH=SINQ*ZCOSHL-COSQ*ZSINHL
      ZN=ZNL
      CC=C1L
      ZE=ZEL
      ZMO=ZMOL
      ASSIGN 40 TO LS
      GO TO 20
   40 SSE = SSE+SE
      SSI=SSI+SI
      SSL=SSL+SL
      SSG=SSG+SGH-COSIQ/SINIQ*SH
      SSH=SSH+SH/SINIQ

*     GEOPOTENTIAL RESONANCE INITIALIZATION FOR 12 HOUR ORBITS

      IRESFL=0
      ISYNFL=0
      IF(XNQ.LT.(.0052359877).AND.XNQ.GT.(.0034906585)) GO TO 70
      IF (XNQ.LT.(8.26E-3) .OR. XNQ.GT.(9.24E-3))    RETURN
      IF (EQ.LT.0.5)    RETURN
      IRESFL =1
      EOC=EQ*EQSQ
      G201=-.306-(EQ-.64)*.440
      IF(EQ.GT.(.65)) GO TO 45
      G211=3.616-13.247*EQ+16.290*EQSQ
      G310=-19.302+117.390*EQ-228.419*EQSQ+156.591*EOC
      G322=-18.9068+109.7927*EQ-214.6334*EQSQ+146.5816*EOC
      G410=-41.122+242.694*EQ-471.094*EQSQ+313.953*EOC
      G422=-146.407+841.880*EQ-1629.014*EQSQ+1083.435*EOC
      G520=-532.114+3017.977*EQ-5740*EQSQ+3708.276*EOC
      GO TO 55
   45 G211=-72.099+331.819*EQ-508.738*EQSQ+266.724*EOC
      G310=-346.844+1582.851*EQ-2415.925*EQSQ+1246.113*EOC
      G322=-342.585+1554.908*EQ-2366.899*EQSQ+1215.972*EOC
      G410=-1052.797+4758.686*EQ-7193.992*EQSQ+3651.957*EOC
      G422=-3581.69+16178.11*EQ-24462.77*EQSQ+12422.52*EOC
      IF(EQ.GT.(.715)) GO TO 50
      G520=1464.74-4664.75*EQ+3763.64*EQSQ
      GO TO 55
   50 G520=-5149.66+29936.92*EQ-54087.36*EQSQ+31324.56*EOC
   55 IF(EQ.GE.(.7)) GO TO 60
      G533=-919.2277+4988.61*EQ-9064.77*EQSQ+5542.21*EOC
      G521 = -822.71072+4568.6173*EQ-8491.4146*EQSQ+5337.524*EOC
      G532 = -853.666+4690.25*EQ-8624.77*EQSQ+5341.4*EOC
      GO TO 65
   60 G533=-37995.78+161616.52*EQ-229838.2*EQSQ+109377.94*EOC
      G521 = -51752.104+218913.95*EQ-309468.16*EQSQ+146349.42*EOC
      G532 = -40023.88+170470.89*EQ-242699.48*EQSQ+115605.82*EOC
   65 SINI2=SINIQ*SINIQ
      F220=.75*(1.+2.*COSIQ+COSQ2)
      F221=1.5*SINI2
      F321=1.875*SINIQ*(1.-2.*COSIQ-3.*COSQ2)
      F322=-1.875*SINIQ*(1.+2.*COSIQ-3.*COSQ2)
      F441=35.*SINI2*F220
      F442=39.3750*SINI2*SINI2
      F522=9.84375*SINIQ*(SINI2*(1.-2.*COSIQ-5.*COSQ2)
     1     +.33333333*(-2.+4.*COSIQ+6.*COSQ2))
      F523 = SINIQ*(4.92187512*SINI2*(-2.-4.*COSIQ+10.*COSQ2)
     *      +6.56250012*(1.+2.*COSIQ-3.*COSQ2))
      F542 = 29.53125*SINIQ*(2.-8.*COSIQ+COSQ2*(-12.+8.*COSIQ
     *      +10.*COSQ2))
      F543=29.53125*SINIQ*(-2.-8.*COSIQ+COSQ2*(12.+8.*COSIQ-10.*COSQ2))
      XNO2=XNQ*XNQ
      AINV2=AQNV*AQNV
      TEMP1 = 3.*XNO2*AINV2
      TEMP = TEMP1*ROOT22
      D2201 = TEMP*F220*G201
      D2211 = TEMP*F221*G211
      TEMP1 = TEMP1*AQNV
      TEMP = TEMP1*ROOT32
      D3210 = TEMP*F321*G310
      D3222 = TEMP*F322*G322
      TEMP1 = TEMP1*AQNV
      TEMP = 2.*TEMP1*ROOT44
      D4410 = TEMP*F441*G410
      D4422 = TEMP*F442*G422
      TEMP1 = TEMP1*AQNV
      TEMP = TEMP1*ROOT52
      D5220 = TEMP*F522*G520
      D5232 = TEMP*F523*G532
      TEMP = 2.*TEMP1*ROOT54
      D5421 = TEMP*F542*G521
      D5433 = TEMP*F543*G533
      XLAMO = XMAO+XNODEO+XNODEO-THGR-THGR
      BFACT = XLLDOT+XNODOT+XNODOT-THDT-THDT
      BFACT=BFACT+SSL+SSH+SSH
      GO TO 80

*      SYNCHRONOUS RESONANCE TERMS INITIALIZATION

   70 IRESFL=1
      ISYNFL=1
      G200=1.0+EQSQ*(-2.5+.8125*EQSQ)
      G310=1.0+2.0*EQSQ
      G300=1.0+EQSQ*(-6.0+6.60937*EQSQ)
      F220=.75*(1.+COSIQ)*(1.+COSIQ)
      F311=.9375*SINIQ*SINIQ*(1.+3.*COSIQ)-.75*(1.+COSIQ)
      F330=1.+COSIQ
      F330=1.875*F330*F330*F330
      DEL1=3.*XNQ*XNQ*AQNV*AQNV
      DEL2=2.*DEL1*F220*G200*Q22
      DEL3=3.*DEL1*F330*G300*Q33*AQNV
      DEL1=DEL1*F311*G310*Q31*AQNV
      FASX2=.13130908
      FASX4=2.8843198
      FASX6=.37448087
      XLAMO=XMAO+XNODEO+OMEGAO-THGR
      BFACT = XLLDOT+XPIDOT-THDT
      BFACT=BFACT+SSL+SSG+SSH
   80 XFACT=BFACT-XNQ
C
C     INITIALIZE INTEGRATOR
C
      XLI=XLAMO
      XNI=XNQ
      ATIME=0.D0
      STEPP=720.D0
      STEPN=-720.D0
      STEP2 = 259200.D0
      RETURN

*     ENTRANCE FOR DEEP SPACE SECULAR EFFECTS

      ENTRY DPSEC(XLL,OMGASM,XNODES,EM,XINC,XN,T)
      XLL=XLL+SSL*T
      OMGASM=OMGASM+SSG*T
      XNODES=XNODES+SSH*T
      EM=EO+SSE*T
      XINC=XINCL+SSI*T
      IF(XINC .GE. 0.) GO TO 90
      XINC = -XINC
      XNODES = XNODES + PI
      OMGASM = OMGASM - PI
   90 IF(IRESFL .EQ. 0) RETURN
  100 IF (ATIME.EQ.0.D0)    GO TO 170
      IF(T.GE.(0.D0).AND.ATIME.LT.(0.D0)) GO TO 170
      IF(T.LT.(0.D0).AND.ATIME.GE.(0.D0)) GO TO 170
  105 IF(DABS(T).GE.DABS(ATIME)) GO TO 120
      DELT=STEPP
      IF (T.GE.0.D0)    DELT = STEPN
  110 ASSIGN 100 TO IRET
      GO TO 160
  120 DELT=STEPN
      IF (T.GT.0.D0)    DELT = STEPP
  125 IF (DABS(T-ATIME).LT.STEPP)    GO TO 130
      ASSIGN 125 TO IRET
      GO TO 160
  130 FT = T-ATIME
      ASSIGN 140 TO IRETN
      GO TO 150
  140 XN = XNI+XNDOT*FT+XNDDT*FT*FT*0.5
      XL = XLI+XLDOT*FT+XNDOT*FT*FT*0.5
      TEMP = -XNODES+THGR+T*THDT
      XLL = XL-OMGASM+TEMP
      IF (ISYNFL.EQ.0)    XLL = XL+TEMP+TEMP
      RETURN
C
C     DOT TERMS CALCULATED
C
  150 IF (ISYNFL.EQ.0)    GO TO 152
      XNDOT=DEL1*SIN (XLI-FASX2)+DEL2*SIN (2.*(XLI-FASX4))
     1     +DEL3*SIN (3.*(XLI-FASX6))
      XNDDT = DEL1*COS(XLI-FASX2)
     *       +2.*DEL2*COS(2.*(XLI-FASX4))
     *       +3.*DEL3*COS(3.*(XLI-FASX6))
      GO TO 154
  152 XOMI = OMEGAQ+OMGDT*ATIME
      X2OMI = XOMI+XOMI
      X2LI = XLI+XLI
      XNDOT = D2201*SIN(X2OMI+XLI-G22)
     *       +D2211*SIN(XLI-G22)
     *       +D3210*SIN(XOMI+XLI-G32)
     *       +D3222*SIN(-XOMI+XLI-G32)
     *       +D4410*SIN(X2OMI+X2LI-G44)
     *       +D4422*SIN(X2LI-G44)
     *       +D5220*SIN(XOMI+XLI-G52)
     *       +D5232*SIN(-XOMI+XLI-G52)
     *       +D5421*SIN(XOMI+X2LI-G54)
     *       +D5433*SIN(-XOMI+X2LI-G54)
      XNDDT = D2201*COS(X2OMI+XLI-G22)
     *       +D2211*COS(XLI-G22)
     *       +D3210*COS(XOMI+XLI-G32)
     *       +D3222*COS(-XOMI+XLI-G32)
     *       +D5220*COS(XOMI+XLI-G52)
     *       +D5232*COS(-XOMI+XLI-G52)
     *       +2.*(D4410*COS(X2OMI+X2LI-G44)
     *       +D4422*COS(X2LI-G44)
     *       +D5421*COS(XOMI+X2LI-G54)
     *       +D5433*COS(-XOMI+X2LI-G54))
  154 XLDOT=XNI+XFACT
      XNDDT = XNDDT*XLDOT
      GO TO IRETN
C
C     INTEGRATOR
C
  160 ASSIGN 165 TO IRETN
      GO TO 150
  165 XLI = XLI+XLDOT*DELT+XNDOT*STEP2
      XNI = XNI+XNDOT*DELT+XNDDT*STEP2
      ATIME=ATIME+DELT
      GO TO IRET
C
C     EPOCH RESTART
C
  170 IF (T.GE.0.D0)    GO TO 175
      DELT=STEPN
      GO TO 180
  175 DELT = STEPP
  180 ATIME = 0.D0
      XNI=XNQ
      XLI=XLAMO
      GO TO 125
C
C     ENTRANCES FOR LUNAR-SOLAR PERIODICS
C
C
      ENTRY DPPER(EM,XINC,OMGASM,XNODES,XLL)
      SINIS = SIN(XINC)
      COSIS = COS(XINC)
      IF (DABS(SAVTSN-T).LT.(30.D0))    GO TO 210
      SAVTSN=T
      ZM=ZMOS+ZNS*T
  205 ZF=ZM+2.*ZES*SIN (ZM)
      SINZF=SIN (ZF)
      F2=.5*SINZF*SINZF-.25
      F3=-.5*SINZF*COS (ZF)
      SES=SE2*F2+SE3*F3
      SIS=SI2*F2+SI3*F3
      SLS=SL2*F2+SL3*F3+SL4*SINZF
      SGHS=SGH2*F2+SGH3*F3+SGH4*SINZF
      SHS=SH2*F2+SH3*F3
      ZM=ZMOL+ZNL*T
      ZF=ZM+2.*ZEL*SIN (ZM)
      SINZF=SIN (ZF)
      F2=.5*SINZF*SINZF-.25
      F3=-.5*SINZF*COS (ZF)
      SEL=EE2*F2+E3*F3
      SIL=XI2*F2+XI3*F3
      SLL=XL2*F2+XL3*F3+XL4*SINZF
      SGHL=XGH2*F2+XGH3*F3+XGH4*SINZF
      SHL=XH2*F2+XH3*F3
      PE=SES+SEL
      PINC=SIS+SIL
      PL=SLS+SLL
  210 PGH=SGHS+SGHL
      PH=SHS+SHL
      XINC = XINC+PINC
      EM = EM+PE
      IF(XQNCL.LT.(.2)) GO TO 220
      GO TO 218
C
C     APPLY PERIODICS DIRECTLY
C
  218 PH=PH/SINIQ
      PGH=PGH-COSIQ*PH
      OMGASM=OMGASM+PGH
      XNODES=XNODES+PH
      XLL = XLL+PL
      GO TO 230
C
C     APPLY PERIODICS WITH LYDDANE MODIFICATION
C
  220 SINOK=SIN(XNODES)
      COSOK=COS(XNODES)
      ALFDP=SINIS*SINOK
      BETDP=SINIS*COSOK
      DALF=PH*COSOK+PINC*COSIS*SINOK
      DBET=-PH*SINOK+PINC*COSIS*COSOK
      ALFDP=ALFDP+DALF
      BETDP=BETDP+DBET
      XLS = XLL+OMGASM+COSIS*XNODES
      DLS=PL+PGH-PINC*XNODES*SINIS
      XLS=XLS+DLS
      XNODES=ACTAN(ALFDP,BETDP)
      XLL = XLL+PL
      OMGASM = XLS-XLL-COS(XINC)*XNODES
  230 CONTINUE
      RETURN
      END
\end{verbatim}
-- DRIVER.FOR --
\begin{verbatim}
*      DRIVER                                                 3 NOV 80

*   WGS-72 PHYSICAL AND GEOPOTENTIAL CONSTANTS
*          CK2= .5*J2*AE**2     CK4=-.375*J4*AE**4

      DOUBLE PRECISION EPOCH,DS50
      COMMON/E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,XNDD6O,BSTAR,
     1            X,Y,Z,XDOT,YDOT,ZDOT,EPOCH,DS50
      COMMON/C1/CK2,CK4,E6A,QOMS2T,S,TOTHRD,
     1           XJ3,XKE,XKMPER,XMNPDA,AE
      COMMON/C2/DE2RA,PI,PIO2,TWOPI,X3PIO2
      DATA IHG/1HG/
      DATA DE2RA,E6A,PI,PIO2,QO,SO,TOTHRD,TWOPI,X3PIO2,XJ2,XJ3,
     1            XJ4,XKE,XKMPER,XMNPDA,AE/.174532925E-1,1.E-6,
     2            3.14159265,1.57079633,120.0,78.0,.66666667,
     4            6.2831853,4.71238898,1.082616E-3,-.253881E-5,
     5            -1.65597E-6,.743669161E-1,6378.135,1440.,1./
      DIMENSION ISET(5)
      CHARACTER ABUF*80(2)
      DATA (ISET(I),I=1,5)/3HSGP,4HSGP4,4HSDP4,4HSGP8,4HSDP8/

*     SELECT EPHEMERIS TYPE AND OUTPUT TIMES

      CK2=.5*XJ2*AE**2
      CK4=-.375*XJ4*AE**4
      QOMS2T=((QO-SO)*AE/XKMPER)**4
      S=AE*(1.+SO/XKMPER)
    2 READ (5,700) IEPT, TS,TF,DELT
      IF(IEPT.LE.0) STOP
      IDEEP=0

*     READ IN MEAN ELEMENTS FROM 2 CARD T(TRANS) OR G(INTERN) FORMAT

      READ (5,706) ABUF
      DECODE(ABUF(1),707)  ITYPE
      IF(ITYPE.EQ.IHG) GO TO 5
      DECODE (ABUF,702) EPOCH,XNDT2O,XNDD6O,IEXP,BSTAR,IBEXP,XINCL,
     1        XNODEO,EO,OMEGAO,XMO,XNO
      GO TO 7
    5 DECODE(ABUF,701) EPOCH,XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,
     1        XNDD6O,IEXP,BSTAR,IBEXP
   7  IF(XNO.LE.0.) STOP
      WRITE(6,704) ABUF,ISET(IEPT)
      IF(IEPT.GT.5) GO TO 900
      XNDD6O=XNDD6O*(10.**IEXP)
      XNODEO=XNODEO*DE2RA
      OMEGAO=OMEGAO*DE2RA
      XMO=XMO*DE2RA
      XINCL=XINCL*DE2RA
      TEMP=TWOPI/XMNPDA/XMNPDA
      XNO=XNO*TEMP*XMNPDA
      XNDT2O=XNDT2O*TEMP
      XNDD6O=XNDD6O*TEMP/XMNPDA

*     INPUT CHECK FOR PERIOD VS EPHEMERIS SELECTED
*     PERIOD GE 225 MINUTES  IS DEEP SPACE

      A1=(XKE/XNO)**TOTHRD
      TEMP=1.5*CK2*(3.*COS(XINCL)**2-1.)/(1.-EO*EO)**1.5
      DEL1=TEMP/(A1*A1)
      AO=A1*(1.-DEL1*(.5*TOTHRD+DEL1*(1.+134./81.*DEL1)))
      DELO=TEMP/(AO*AO)
      XNODP=XNO/(1.+DELO)
      IF((TWOPI/XNODP/XMNPDA) .GE. .15625) IDEEP=1

      BSTAR=BSTAR*(10.**IBEXP)/AE
      TSINCE=TS
      IFLAG=1
      IF(IDEEP .EQ. 1 .AND. (IEPT .EQ. 1 .OR. IEPT .EQ. 2
     1           .OR. IEPT .EQ. 4)) GO TO 800
    9 IF(IDEEP .EQ. 0 .AND. (IEPT .EQ. 3 .OR. IEPT .EQ. 5))
     1           GO TO 850
   10 GO TO (21,22,23,24,25), IEPT
   21 CALL SGP(IFLAG,TSINCE)
      GO TO 60
   22 CALL SGP4(IFLAG,TSINCE)
      GO TO 60
   23 CALL SDP4(IFLAG,TSINCE)
      GO TO 60
   24 CALL SGP8(IFLAG,TSINCE)
      GO TO 60
   25 CALL SDP8(IFLAG,TSINCE)
   60 X=X*XKMPER/AE
      Y=Y*XKMPER/AE
      Z=Z*XKMPER/AE
      XDOT=XDOT*XKMPER/AE*XMNPDA/86400.
      YDOT=YDOT*XKMPER/AE*XMNPDA/86400.
      ZDOT=ZDOT*XKMPER/AE*XMNPDA/86400.
      WRITE(6,705) TSINCE,X,Y,Z,XDOT,YDOT,ZDOT
      TSINCE=TSINCE+DELT
      IF(ABS(TSINCE) .GT. ABS(TF)) GO TO 2
      GO TO 10
  700 FORMAT(I1,3F10.0)
  701 FORMAT(29X,D14.8,1X,3F8.4,/,6X,F8.7,F8.4,1X,2F11.9,1X,F6.5,I2,
     1       4X,F8.7,I2)
  702 FORMAT(18X,D14.8,1X,F10.8,2(1X,F6.5,I2),/,7X,2(1X,F8.4),1X,
     1       F7.7,2(1X,F8.4),1X,F11.8)
  703 FORMAT(79X,A1)
  704 FORMAT(1H1,A80,/,1X,A80,//,1X,A4,7H TSINCE,
     1 14X,1HX,16X,1HY,16X,1HZ,14X,
     1 4HXDOT,13X,4HYDOT,13X,4HZDOT,//)
  705 FORMAT(7F17.8)
  706 FORMAT(A80)
  707 FORMAT(79X,A1)
  930 FORMAT("SHOULD USE DEEP SPACE EPHEMERIS")
  940 FORMAT("SHOULD USE NEAR EARTH EPHEMERIS")
  950 FORMAT("EPHEMERIS NUMBER",I2," NOT LEGAL, WILL SKIP THIS CASE")
  800 WRITE(6,930)
      GO TO 9
  850 WRITE(6,940)
      GO TO 10
  900 WRITE(6,950) IEPT
      GO TO 2
      END
\end{verbatim}
-- FMOD2P.FOR --
\begin{verbatim}
      FUNCTION FMOD2P(X)
      COMMON/C2/DE2RA,PI,PIO2,TWOPI,X3PIO2
      FMOD2P=X
      I=FMOD2P/TWOPI
       FMOD2P=FMOD2P-I*TWOPI
      IF(FMOD2P.LT.0) FMOD2P=FMOD2P+TWOPI
      RETURN
      END
\end{verbatim}
-- SGP.FOR --
\begin{verbatim}
*      SGP                                             31 OCT 80
      SUBROUTINE SGP(IFLAG,TSINCE)
      COMMON/E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,XNDD6O,BSTAR,
     1         X,Y,Z,XDOT,YDOT,ZDOT,EPOCH,DS50
      COMMON/C1/CK2,CK4,E6A,QOMS2T,S,TOTHRD,
     1        XJ3,XKE,XKMPER,XMNPDA,AE
      DOUBLE PRECISION EPOCH, DS50

      IF(IFLAG.EQ.0) GO TO 19

*      INITIALIZATION

      C1= CK2*1.5
      C2= CK2/4.0
      C3= CK2/2.0
      C4= XJ3*AE**3/(4.0*CK2)
      COSIO=COS(XINCL)
      SINIO=SIN(XINCL)
      A1=(XKE/XNO)**TOTHRD
      D1=     C1/A1/A1*(3.*COSIO*COSIO-1.)/(1.-EO*EO)**1.5
      AO=A1*(1.-1./3.*D1-D1*D1-134./81.*D1*D1*D1)
      PO=AO*(1.-EO*EO)
      QO=AO*(1.-EO)
      XLO=XMO+OMEGAO+XNODEO
      D1O= C3 *SINIO*SINIO
      D2O= C2 *(7.*COSIO*COSIO-1.)
      D3O=C1*COSIO
      D4O=D3O*SINIO
      PO2NO=XNO/(PO*PO)
      OMGDT=C1*PO2NO*(5.*COSIO*COSIO-1.)
      XNODOT=-2.*D3O*PO2NO
      C5=.5*C4*SINIO*(3.+5.*COSIO)/(1.+COSIO)
      C6=C4*SINIO
      IFLAG=0

*      UPDATE FOR SECULAR GRAVITY AND ATMOSPHERIC DRAG

   19 A=XNO+(2.*XNDT2O+3.*XNDD6O*TSINCE)*TSINCE
      A=AO*(XNO/A)**TOTHRD
      E=E6A
      IF(A.GT.QO) E=1.-QO/A
      P=A*(1.-E*E)
      XNODES= XNODEO+XNODOT*TSINCE
      OMGAS= OMEGAO+OMGDT*TSINCE
      XLS=FMOD2P(XLO+(XNO+OMGDT+XNODOT+(XNDT2O+XNDD6O*TSINCE)*
     1 TSINCE)*TSINCE)

*      LONG PERIOD PERIODICS

      AXNSL=E*COS(OMGAS)
      AYNSL=E*SIN(OMGAS)-C6/P
      XL=FMOD2P(XLS-C5/P*AXNSL)

*      SOLVE KEPLERS EQUATION

      U=FMOD2P(XL-XNODES)
      ITEM3=0
      EO1=U
      TEM5=1.
   20 SINEO1=SIN(EO1)
      COSEO1=COS(EO1)
      IF(ABS(TEM5).LT.E6A) GO TO 30
      IF(ITEM3.GE.10) GO TO 30
      ITEM3=ITEM3+1
      TEM5=1.-COSEO1*AXNSL-SINEO1*AYNSL
      TEM5=(U-AYNSL*COSEO1+AXNSL*SINEO1-EO1)/TEM5
      TEM2=ABS(TEM5)
      IF(TEM2.GT.1.) TEM5=TEM2/TEM5
      EO1=EO1+TEM5
      GO TO 20

*      SHORT PERIOD PRELIMINARY QUANTITIES

   30 ECOSE=AXNSL*COSEO1+AYNSL*SINEO1
      ESINE=AXNSL*SINEO1-AYNSL*COSEO1
      EL2=AXNSL*AXNSL+AYNSL*AYNSL
      PL=A*(1.-EL2)
      PL2=PL*PL
      R=A*(1.-ECOSE)
      RDOT=XKE*SQRT(A)/R*ESINE
      RVDOT=XKE*SQRT(PL)/R
      TEMP=ESINE/(1.+SQRT(1.-EL2))
      SINU=A/R*(SINEO1-AYNSL-AXNSL*TEMP)
      COSU=A/R*(COSEO1-AXNSL+AYNSL*TEMP)
      SU=ACTAN(SINU,COSU)

*      UPDATE FOR SHORT PERIODICS

      SIN2U=(COSU+COSU)*SINU
      COS2U=1.-2.*SINU*SINU
      RK=R+D1O/PL*COS2U
      UK=SU-D2O/PL2*SIN2U
      XNODEK=XNODES+D3O*SIN2U/PL2
      XINCK =XINCL+D4O/PL2*COS2U

*      ORIENTATION VECTORS

      SINUK=SIN(UK)
      COSUK=COS(UK)
      SINNOK=SIN(XNODEK)
      COSNOK=COS(XNODEK)
      SINIK=SIN(XINCK)
      COSIK=COS(XINCK)
      XMX=-SINNOK*COSIK
      XMY=COSNOK*COSIK
      UX=XMX*SINUK+COSNOK*COSUK
      UY=XMY*SINUK+SINNOK*COSUK
      UZ=SINIK*SINUK
      VX=XMX*COSUK-COSNOK*SINUK
      VY=XMY*COSUK-SINNOK*SINUK
      VZ=SINIK*COSUK

*      POSITION AND VELOCITY

      X=RK*UX
      Y=RK*UY
      Z=RK*UZ
      XDOT=RDOT*UX
      YDOT=RDOT*UY
      ZDOT=RDOT*UZ
      XDOT=RVDOT*VX+XDOT
      YDOT=RVDOT*VY+YDOT
      ZDOT=RVDOT*VZ+ZDOT

      RETURN
      END
\end{verbatim}
-- SGP4.FOR --
\begin{verbatim}
*      SGP4                                               3 NOV 80
      SUBROUTINE SGP4(IFLAG,TSINCE)
      COMMON/E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,
     1           XNDD6O,BSTAR,X,Y,Z,XDOT,YDOT,ZDOT,EPOCH,DS50
      COMMON/C1/CK2,CK4,E6A,QOMS2T,S,TOTHRD,
     1           XJ3,XKE,XKMPER,XMNPDA,AE
      DOUBLE PRECISION EPOCH, DS50

      IF (IFLAG .EQ. 0) GO TO 100

*      RECOVER ORIGINAL MEAN MOTION (XNODP) AND SEMIMAJOR AXIS (AODP)
*      FROM INPUT ELEMENTS

      A1=(XKE/XNO)**TOTHRD
      COSIO=COS(XINCL)
      THETA2=COSIO*COSIO
      X3THM1=3.*THETA2-1.
      EOSQ=EO*EO
      BETAO2=1.-EOSQ
      BETAO=SQRT(BETAO2)
      DEL1=1.5*CK2*X3THM1/(A1*A1*BETAO*BETAO2)
      AO=A1*(1.-DEL1*(.5*TOTHRD+DEL1*(1.+134./81.*DEL1)))
      DELO=1.5*CK2*X3THM1/(AO*AO*BETAO*BETAO2)
      XNODP=XNO/(1.+DELO)
      AODP=AO/(1.-DELO)

*      INITIALIZATION

*      FOR PERIGEE LESS THAN 220 KILOMETERS, THE ISIMP FLAG IS SET AND
*      THE EQUATIONS ARE TRUNCATED TO LINEAR VARIATION IN SQRT A AND
*      QUADRATIC VARIATION IN MEAN ANOMALY.  ALSO, THE C3 TERM, THE
*      DELTA OMEGA TERM, AND THE DELTA M TERM ARE DROPPED.

      ISIMP=0
      IF((AODP*(1.-EO)/AE) .LT. (220./XKMPER+AE)) ISIMP=1

*      FOR PERIGEE BELOW 156 KM, THE VALUES OF
*      S AND QOMS2T ARE ALTERED

      S4=S
      QOMS24=QOMS2T
      PERIGE=(AODP*(1.-EO)-AE)*XKMPER
      IF(PERIGE .GE. 156.) GO TO 10
      S4=PERIGE-78.
      IF(PERIGE .GT. 98.) GO TO 9
      S4=20.
    9 QOMS24=((120.-S4)*AE/XKMPER)**4
      S4=S4/XKMPER+AE
   10 PINVSQ=1./(AODP*AODP*BETAO2*BETAO2)
      TSI=1./(AODP-S4)
      ETA=AODP*EO*TSI
      ETASQ=ETA*ETA
      EETA=EO*ETA
      PSISQ=ABS(1.-ETASQ)
      COEF=QOMS24*TSI**4
      COEF1=COEF/PSISQ**3.5
      C2=COEF1*XNODP*(AODP*(1.+1.5*ETASQ+EETA*(4.+ETASQ))+.75*
     1         CK2*TSI/PSISQ*X3THM1*(8.+3.*ETASQ*(8.+ETASQ)))
      C1=BSTAR*C2
      SINIO=SIN(XINCL)
      A3OVK2=-XJ3/CK2*AE**3
      C3=COEF*TSI*A3OVK2*XNODP*AE*SINIO/EO
      X1MTH2=1.-THETA2
      C4=2.*XNODP*COEF1*AODP*BETAO2*(ETA*
     1         (2.+.5*ETASQ)+EO*(.5+2.*ETASQ)-2.*CK2*TSI/
     2         (AODP*PSISQ)*(-3.*X3THM1*(1.-2.*EETA+ETASQ*
     3         (1.5-.5*EETA))+.75*X1MTH2*(2.*ETASQ-EETA*
     4         (1.+ETASQ))*COS(2.*OMEGAO)))
      C5=2.*COEF1*AODP*BETAO2*(1.+2.75*(ETASQ+EETA)+EETA*ETASQ)
      THETA4=THETA2*THETA2
      TEMP1=3.*CK2*PINVSQ*XNODP
      TEMP2=TEMP1*CK2*PINVSQ
      TEMP3=1.25*CK4*PINVSQ*PINVSQ*XNODP
      XMDOT=XNODP+.5*TEMP1*BETAO*X3THM1+.0625*TEMP2*BETAO*
     1         (13.-78.*THETA2+137.*THETA4)
      X1M5TH=1.-5.*THETA2
      OMGDOT=-.5*TEMP1*X1M5TH+.0625*TEMP2*(7.-114.*THETA2+
     1         395.*THETA4)+TEMP3*(3.-36.*THETA2+49.*THETA4)
      XHDOT1=-TEMP1*COSIO
      XNODOT=XHDOT1+(.5*TEMP2*(4.-19.*THETA2)+2.*TEMP3*(3.-
     1         7.*THETA2))*COSIO
      OMGCOF=BSTAR*C3*COS(OMEGAO)
      XMCOF=-TOTHRD*COEF*BSTAR*AE/EETA
      XNODCF=3.5*BETAO2*XHDOT1*C1
      T2COF=1.5*C1
      XLCOF=.125*A3OVK2*SINIO*(3.+5.*COSIO)/(1.+COSIO)
      AYCOF=.25*A3OVK2*SINIO
      DELMO=(1.+ETA*COS(XMO))**3
      SINMO=SIN(XMO)
      X7THM1=7.*THETA2-1.
      IF(ISIMP .EQ. 1) GO TO 90
      C1SQ=C1*C1
      D2=4.*AODP*TSI*C1SQ
      TEMP=D2*TSI*C1/3.
      D3=(17.*AODP+S4)*TEMP
      D4=.5*TEMP*AODP*TSI*(221.*AODP+31.*S4)*C1
      T3COF=D2+2.*C1SQ
      T4COF=.25*(3.*D3+C1*(12.*D2+10.*C1SQ))
      T5COF=.2*(3.*D4+12.*C1*D3+6.*D2*D2+15.*C1SQ*(
     1         2.*D2+C1SQ))
   90 IFLAG=0

*      UPDATE FOR SECULAR GRAVITY AND ATMOSPHERIC DRAG

  100 XMDF=XMO+XMDOT*TSINCE
      OMGADF=OMEGAO+OMGDOT*TSINCE
      XNODDF=XNODEO+XNODOT*TSINCE
      OMEGA=OMGADF
      XMP=XMDF
      TSQ=TSINCE*TSINCE
      XNODE=XNODDF+XNODCF*TSQ
      TEMPA=1.-C1*TSINCE
      TEMPE=BSTAR*C4*TSINCE
      TEMPL=T2COF*TSQ
      IF(ISIMP .EQ. 1) GO TO 110
      DELOMG=OMGCOF*TSINCE
      DELM=XMCOF*((1.+ETA*COS(XMDF))**3-DELMO)
      TEMP=DELOMG+DELM
      XMP=XMDF+TEMP
      OMEGA=OMGADF-TEMP
      TCUBE=TSQ*TSINCE
      TFOUR=TSINCE*TCUBE
      TEMPA=TEMPA-D2*TSQ-D3*TCUBE-D4*TFOUR
      TEMPE=TEMPE+BSTAR*C5*(SIN(XMP)-SINMO)
      TEMPL=TEMPL+T3COF*TCUBE+
     1         TFOUR*(T4COF+TSINCE*T5COF)
  110 A=AODP*TEMPA**2
      E=EO-TEMPE
      XL=XMP+OMEGA+XNODE+XNODP*TEMPL
      BETA=SQRT(1.-E*E)
      XN=XKE/A**1.5

*      LONG PERIOD PERIODICS

      AXN=E*COS(OMEGA)
      TEMP=1./(A*BETA*BETA)
      XLL=TEMP*XLCOF*AXN
      AYNL=TEMP*AYCOF
      XLT=XL+XLL
      AYN=E*SIN(OMEGA)+AYNL

*      SOLVE KEPLERS EQUATION

      CAPU=FMOD2P(XLT-XNODE)
      TEMP2=CAPU
      DO 130 I=1,10
      SINEPW=SIN(TEMP2)
      COSEPW=COS(TEMP2)
      TEMP3=AXN*SINEPW
      TEMP4=AYN*COSEPW
      TEMP5=AXN*COSEPW
      TEMP6=AYN*SINEPW
      EPW=(CAPU-TEMP4+TEMP3-TEMP2)/(1.-TEMP5-TEMP6)+TEMP2
      IF(ABS(EPW-TEMP2) .LE. E6A) GO TO 140
  130 TEMP2=EPW

*      SHORT PERIOD PRELIMINARY QUANTITIES

  140 ECOSE=TEMP5+TEMP6
      ESINE=TEMP3-TEMP4
      ELSQ=AXN*AXN+AYN*AYN
      TEMP=1.-ELSQ
      PL=A*TEMP
      R=A*(1.-ECOSE)
      TEMP1=1./R
      RDOT=XKE*SQRT(A)*ESINE*TEMP1
      RFDOT=XKE*SQRT(PL)*TEMP1
      TEMP2=A*TEMP1
      BETAL=SQRT(TEMP)
      TEMP3=1./(1.+BETAL)
      COSU=TEMP2*(COSEPW-AXN+AYN*ESINE*TEMP3)
      SINU=TEMP2*(SINEPW-AYN-AXN*ESINE*TEMP3)
      U=ACTAN(SINU,COSU)
      SIN2U=2.*SINU*COSU
      COS2U=2.*COSU*COSU-1.
      TEMP=1./PL
      TEMP1=CK2*TEMP
      TEMP2=TEMP1*TEMP

*      UPDATE FOR SHORT PERIODICS

      RK=R*(1.-1.5*TEMP2*BETAL*X3THM1)+.5*TEMP1*X1MTH2*COS2U
      UK=U-.25*TEMP2*X7THM1*SIN2U
      XNODEK=XNODE+1.5*TEMP2*COSIO*SIN2U
      XINCK=XINCL+1.5*TEMP2*COSIO*SINIO*COS2U
      RDOTK=RDOT-XN*TEMP1*X1MTH2*SIN2U
      RFDOTK=RFDOT+XN*TEMP1*(X1MTH2*COS2U+1.5*X3THM1)

*      ORIENTATION VECTORS

      SINUK=SIN(UK)
      COSUK=COS(UK)
      SINIK=SIN(XINCK)
      COSIK=COS(XINCK)
      SINNOK=SIN(XNODEK)
      COSNOK=COS(XNODEK)
      XMX=-SINNOK*COSIK
      XMY=COSNOK*COSIK
      UX=XMX*SINUK+COSNOK*COSUK
      UY=XMY*SINUK+SINNOK*COSUK
      UZ=SINIK*SINUK
      VX=XMX*COSUK-COSNOK*SINUK
      VY=XMY*COSUK-SINNOK*SINUK
      VZ=SINIK*COSUK

*      POSITION AND VELOCITY

      X=RK*UX
      Y=RK*UY
      Z=RK*UZ
      XDOT=RDOTK*UX+RFDOTK*VX
      YDOT=RDOTK*UY+RFDOTK*VY
      ZDOT=RDOTK*UZ+RFDOTK*VZ

      RETURN
      END
\end{verbatim}
-- SDP4.FOR --
\begin{verbatim}
*      SDP4                                               3 NOV 80
      SUBROUTINE SDP4(IFLAG,TSINCE)
      COMMON/E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,
     1           XNDD6O,BSTAR,X,Y,Z,XDOT,YDOT,ZDOT,EPOCH,DS50
      COMMON/C1/CK2,CK4,E6A,QOMS2T,S,TOTHRD,
     1           XJ3,XKE,XKMPER,XMNPDA,AE
      DOUBLE PRECISION EPOCH, DS50

      IF  (IFLAG .EQ. 0) GO TO 100

*      RECOVER ORIGINAL MEAN MOTION (XNODP) AND SEMIMAJOR AXIS (AODP)
*      FROM INPUT ELEMENTS

      A1=(XKE/XNO)**TOTHRD
      COSIO=COS(XINCL)
      THETA2=COSIO*COSIO
      X3THM1=3.*THETA2-1.
      EOSQ=EO*EO
      BETAO2=1.-EOSQ
      BETAO=SQRT(BETAO2)
      DEL1=1.5*CK2*X3THM1/(A1*A1*BETAO*BETAO2)
      AO=A1*(1.-DEL1*(.5*TOTHRD+DEL1*(1.+134./81.*DEL1)))
      DELO=1.5*CK2*X3THM1/(AO*AO*BETAO*BETAO2)
      XNODP=XNO/(1.+DELO)
      AODP=AO/(1.-DELO)

*      INITIALIZATION

*      FOR PERIGEE BELOW 156 KM, THE VALUES OF
*      S AND QOMS2T ARE ALTERED

      S4=S
      QOMS24=QOMS2T
      PERIGE=(AODP*(1.-EO)-AE)*XKMPER
      IF(PERIGE .GE. 156.) GO TO 10
      S4=PERIGE-78.
      IF(PERIGE .GT. 98.) GO TO 9
      S4=20.
    9 QOMS24=((120.-S4)*AE/XKMPER)**4
      S4=S4/XKMPER+AE
   10 PINVSQ=1./(AODP*AODP*BETAO2*BETAO2)
      SING=SIN(OMEGAO)
      COSG=COS(OMEGAO)
      TSI=1./(AODP-S4)
      ETA=AODP*EO*TSI
      ETASQ=ETA*ETA
      EETA=EO*ETA
      PSISQ=ABS(1.-ETASQ)
      COEF=QOMS24*TSI**4
      COEF1=COEF/PSISQ**3.5
      C2=COEF1*XNODP*(AODP*(1.+1.5*ETASQ+EETA*(4.+ETASQ))+.75*
     1         CK2*TSI/PSISQ*X3THM1*(8.+3.*ETASQ*(8.+ETASQ)))
      C1=BSTAR*C2
      SINIO=SIN(XINCL)
      A3OVK2=-XJ3/CK2*AE**3
      X1MTH2=1.-THETA2
      C4=2.*XNODP*COEF1*AODP*BETAO2*(ETA*
     1         (2.+.5*ETASQ)+EO*(.5+2.*ETASQ)-2.*CK2*TSI/
     2         (AODP*PSISQ)*(-3.*X3THM1*(1.-2.*EETA+ETASQ*
     3         (1.5-.5*EETA))+.75*X1MTH2*(2.*ETASQ-EETA*
     4         (1.+ETASQ))*COS(2.*OMEGAO)))
      THETA4=THETA2*THETA2
      TEMP1=3.*CK2*PINVSQ*XNODP
      TEMP2=TEMP1*CK2*PINVSQ
      TEMP3=1.25*CK4*PINVSQ*PINVSQ*XNODP
      XMDOT=XNODP+.5*TEMP1*BETAO*X3THM1+.0625*TEMP2*BETAO*
     1         (13.-78.*THETA2+137.*THETA4)
      X1M5TH=1.-5.*THETA2
      OMGDOT=-.5*TEMP1*X1M5TH+.0625*TEMP2*(7.-114.*THETA2+
     1         395.*THETA4)+TEMP3*(3.-36.*THETA2+49.*THETA4)
      XHDOT1=-TEMP1*COSIO
      XNODOT=XHDOT1+(.5*TEMP2*(4.-19.*THETA2)+2.*TEMP3*(3.-
     1         7.*THETA2))*COSIO
      XNODCF=3.5*BETAO2*XHDOT1*C1
      T2COF=1.5*C1
      XLCOF=.125*A3OVK2*SINIO*(3.+5.*COSIO)/(1.+COSIO)
      AYCOF=.25*A3OVK2*SINIO
      X7THM1=7.*THETA2-1.
   90 IFLAG=0
      CALL DPINIT(EOSQ,SINIO,COSIO,BETAO,AODP,THETA2,
     1         SING,COSG,BETAO2,XMDOT,OMGDOT,XNODOT,XNODP)

*      UPDATE FOR SECULAR GRAVITY AND ATMOSPHERIC DRAG

  100 XMDF=XMO+XMDOT*TSINCE
      OMGADF=OMEGAO+OMGDOT*TSINCE
      XNODDF=XNODEO+XNODOT*TSINCE
      TSQ=TSINCE*TSINCE
      XNODE=XNODDF+XNODCF*TSQ
      TEMPA=1.-C1*TSINCE
      TEMPE=BSTAR*C4*TSINCE
      TEMPL=T2COF*TSQ
      XN=XNODP
      CALL DPSEC(XMDF,OMGADF,XNODE,EM,XINC,XN,TSINCE)
      A=(XKE/XN)**TOTHRD*TEMPA**2
      E=EM-TEMPE
      XMAM=XMDF+XNODP*TEMPL
      CALL DPPER(E,XINC,OMGADF,XNODE,XMAM)
      XL=XMAM+OMGADF+XNODE
      BETA=SQRT(1.-E*E)
      XN=XKE/A**1.5

*      LONG PERIOD PERIODICS

      AXN=E*COS(OMGADF)
      TEMP=1./(A*BETA*BETA)
      XLL=TEMP*XLCOF*AXN
      AYNL=TEMP*AYCOF
      XLT=XL+XLL
      AYN=E*SIN(OMGADF)+AYNL

*      SOLVE KEPLERS EQUATION

      CAPU=FMOD2P(XLT-XNODE)
      TEMP2=CAPU
      DO 130 I=1,10
      SINEPW=SIN(TEMP2)
      COSEPW=COS(TEMP2)
      TEMP3=AXN*SINEPW
      TEMP4=AYN*COSEPW
      TEMP5=AXN*COSEPW
      TEMP6=AYN*SINEPW
      EPW=(CAPU-TEMP4+TEMP3-TEMP2)/(1.-TEMP5-TEMP6)+TEMP2
      IF(ABS(EPW-TEMP2) .LE. E6A) GO TO 140
  130 TEMP2=EPW

*      SHORT PERIOD PRELIMINARY QUANTITIES

  140 ECOSE=TEMP5+TEMP6
      ESINE=TEMP3-TEMP4
      ELSQ=AXN*AXN+AYN*AYN
      TEMP=1.-ELSQ
      PL=A*TEMP
      R=A*(1.-ECOSE)
      TEMP1=1./R
      RDOT=XKE*SQRT(A)*ESINE*TEMP1
      RFDOT=XKE*SQRT(PL)*TEMP1
      TEMP2=A*TEMP1
      BETAL=SQRT(TEMP)
      TEMP3=1./(1.+BETAL)
      COSU=TEMP2*(COSEPW-AXN+AYN*ESINE*TEMP3)
      SINU=TEMP2*(SINEPW-AYN-AXN*ESINE*TEMP3)
      U=ACTAN(SINU,COSU)
      SIN2U=2.*SINU*COSU
      COS2U=2.*COSU*COSU-1.
      TEMP=1./PL
      TEMP1=CK2*TEMP
      TEMP2=TEMP1*TEMP

*      UPDATE FOR SHORT PERIODICS

      RK=R*(1.-1.5*TEMP2*BETAL*X3THM1)+.5*TEMP1*X1MTH2*COS2U
      UK=U-.25*TEMP2*X7THM1*SIN2U
      XNODEK=XNODE+1.5*TEMP2*COSIO*SIN2U
      XINCK=XINC+1.5*TEMP2*COSIO*SINIO*COS2U
      RDOTK=RDOT-XN*TEMP1*X1MTH2*SIN2U
      RFDOTK=RFDOT+XN*TEMP1*(X1MTH2*COS2U+1.5*X3THM1)

*      ORIENTATION VECTORS

      SINUK=SIN(UK)
      COSUK=COS(UK)
      SINIK=SIN(XINCK)
      COSIK=COS(XINCK)
      SINNOK=SIN(XNODEK)
      COSNOK=COS(XNODEK)
      XMX=-SINNOK*COSIK
      XMY=COSNOK*COSIK
      UX=XMX*SINUK+COSNOK*COSUK
      UY=XMY*SINUK+SINNOK*COSUK
      UZ=SINIK*SINUK
      VX=XMX*COSUK-COSNOK*SINUK
      VY=XMY*COSUK-SINNOK*SINUK
      VZ=SINIK*COSUK

*      POSITION AND VELOCITY

      X=RK*UX
      Y=RK*UY
      Z=RK*UZ
      XDOT=RDOTK*UX+RFDOTK*VX
      YDOT=RDOTK*UY+RFDOTK*VY
      ZDOT=RDOTK*UZ+RFDOTK*VZ

      RETURN
      END
\end{verbatim}
-- SGP8.FOR --
\begin{verbatim}
*      SGP8                                              14 NOV 80
      SUBROUTINE SGP8(IFLAG,TSINCE)
      COMMON/E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,
     1           XNDD6O,BSTAR,X,Y,Z,XDOT,YDOT,ZDOT,EPOCH,DS50
      COMMON/C1/CK2,CK4,E6A,QOMS2T,S,TOTHRD,
     1           XJ3,XKE,XKMPER,XMNPDA,AE
      DOUBLE PRECISION EPOCH, DS50
      DATA RHO/.15696615/

      IF (IFLAG .EQ. 0) GO TO 100

*      RECOVER ORIGINAL MEAN MOTION (XNODP) AND SEMIMAJOR AXIS (AODP)
*      FROM INPUT ELEMENTS --------- CALCULATE BALLISTIC COEFFICIENT
*      (B TERM) FROM INPUT B* DRAG TERM

      A1=(XKE/XNO)**TOTHRD
      COSI=COS(XINCL)
      THETA2=COSI*COSI
      TTHMUN=3.*THETA2-1.
      EOSQ=EO*EO
      BETAO2=1.-EOSQ
      BETAO=SQRT(BETAO2)
      DEL1=1.5*CK2*TTHMUN/(A1*A1*BETAO*BETAO2)
      AO=A1*(1.-DEL1*(.5*TOTHRD+DEL1*(1.+134./81.*DEL1)))
      DELO=1.5*CK2*TTHMUN/(AO*AO*BETAO*BETAO2)
      AODP=AO/(1.-DELO)
      XNODP=XNO/(1.+DELO)
      B=2.*BSTAR/RHO

*      INITIALIZATION

      ISIMP=0
      PO=AODP*BETAO2
      POM2=1./(PO*PO)
      SINI=SIN(XINCL)
      SING=SIN(OMEGAO)
      COSG=COS(OMEGAO)
      TEMP=.5*XINCL
      SINIO2=SIN(TEMP)
      COSIO2=COS(TEMP)
      THETA4=THETA2**2
      UNM5TH=1.-5.*THETA2
      UNMTH2=1.-THETA2
      A3COF=-XJ3/CK2*AE**3
      PARDT1=3.*CK2*POM2*XNODP
      PARDT2=PARDT1*CK2*POM2
      PARDT4=1.25*CK4*POM2*POM2*XNODP
      XMDT1=.5*PARDT1*BETAO*TTHMUN
      XGDT1=-.5*PARDT1*UNM5TH
      XHDT1=-PARDT1*COSI
      XLLDOT=XNODP+XMDT1+
     2           .0625*PARDT2*BETAO*(13.-78.*THETA2+137.*THETA4)
      OMGDT=XGDT1+
     1      .0625*PARDT2*(7.-114.*THETA2+395.*THETA4)+PARDT4*(3.-36.*
     2         THETA2+49.*THETA4)
      XNODOT=XHDT1+
     1       (.5*PARDT2*(4.-19.*THETA2)+2.*PARDT4*(3.-7.*THETA2))*COSI
      TSI=1./(PO-S)
      ETA=EO*S*TSI
      ETA2=ETA**2
      PSIM2=ABS(1./(1.-ETA2))
      ALPHA2=1.+EOSQ
      EETA=EO*ETA
      COS2G=2.*COSG**2-1.
      D5=TSI*PSIM2
      D1=D5/PO
      D2=12.+ETA2*(36.+4.5*ETA2)
      D3=ETA2*(15.+2.5*ETA2)
      D4=ETA*(5.+3.75*ETA2)
      B1=CK2*TTHMUN
      B2=-CK2*UNMTH2
      B3=A3COF*SINI
      C0=.5*B*RHO*QOMS2T*XNODP*AODP*TSI**4*PSIM2**3.5/SQRT(ALPHA2)
      C1=1.5*XNODP*ALPHA2**2*C0
      C4=D1*D3*B2
      C5=D5*D4*B3
      XNDT=C1*(
     1  (2.+ETA2*(3.+34.*EOSQ)+5.*EETA*(4.+ETA2)+8.5*EOSQ)+
     1   D1*D2*B1+   C4*COS2G+C5*SING)
      XNDTN=XNDT/XNODP

*      IF DRAG IS VERY SMALL, THE ISIMP FLAG IS SET AND THE
*      EQUATIONS ARE TRUNCATED TO LINEAR VARIATION IN MEAN
*      MOTION AND QUADRATIC VARIATION IN MEAN ANOMALY

      IF(ABS(XNDTN*XMNPDA) .LT. 2.16E-3) GO TO 50
      D6=ETA*(30.+22.5*ETA2)
      D7=ETA*(5.+12.5*ETA2)
      D8=1.+ETA2*(6.75+ETA2)
      C8=D1*D7*B2
      C9=D5*D8*B3
      EDOT=-C0*(
     1   ETA*(4.+ETA2+EOSQ*(15.5+7.*ETA2))+EO*(5.+15.*ETA2)+
     1   D1*D6*B1 +
     1   C8*COS2G+C9*SING)
      D20=.5*TOTHRD*XNDTN
      ALDTAL=EO*EDOT/ALPHA2
      TSDTTS=2.*AODP*TSI*(D20*BETAO2+EO*EDOT)
      ETDT=(EDOT+EO*TSDTTS)*TSI*S
      PSDTPS=-ETA*ETDT*PSIM2
      SIN2G=2.*SING*COSG
      C0DTC0=D20+4.*TSDTTS-ALDTAL-7.*PSDTPS
      C1DTC1=XNDTN+4.*ALDTAL+C0DTC0
      D9=ETA*(6.+68.*EOSQ)+EO*(20.+15.*ETA2)
      D10=5.*ETA*(4.+ETA2)+EO*(17.+68.*ETA2)
      D11=ETA*(72.+18.*ETA2)
      D12=ETA*(30.+10.*ETA2)
      D13=5.+11.25*ETA2
      D14=TSDTTS-2.*PSDTPS
      D15=2.*(D20+EO*EDOT/BETAO2)
      D1DT=D1*(D14+D15)
      D2DT=ETDT*D11
      D3DT=ETDT*D12
      D4DT=ETDT*D13
      D5DT=D5*D14
      C4DT=B2*(D1DT*D3+D1*D3DT)
      C5DT=B3*(D5DT*D4+D5*D4DT)
      D16=
     1     D9*ETDT+D10*EDOT +
     1     B1*(D1DT*D2+D1*D2DT) +
     1     C4DT*COS2G+C5DT*SING+XGDT1*(C5*COSG-2.*C4*SIN2G)
      XNDDT=C1DTC1*XNDT+C1*D16
      EDDOT=C0DTC0*EDOT-C0*(
     1     (4.+3.*ETA2+30.*EETA+EOSQ*(15.5+21.*ETA2))*ETDT+(5.+15.*ETA2
     '         +EETA*(31.+14.*ETA2))*EDOT +
     1     B1*(D1DT*D6+D1*ETDT*(30.+67.5*ETA2))  +
     1     B2*(D1DT*D7+D1*ETDT*(5.+37.5*ETA2))*COS2G+
     1     B3*(D5DT*D8+D5*ETDT*ETA*(13.5+4.*ETA2))*SING+XGDT1*(C9*
     '         COSG-2.*C8*SIN2G))
      D25=EDOT**2
      D17=XNDDT/XNODP-XNDTN**2
      TSDDTS=2.*TSDTTS*(TSDTTS-D20)+AODP*TSI*(TOTHRD*BETAO2*D17-4.*D20*
     '         EO*EDOT+2.*(D25+EO*EDDOT))
      ETDDT =(EDDOT+2.*EDOT*TSDTTS)*TSI*S+TSDDTS*ETA
      D18=TSDDTS-TSDTTS**2
      D19=-PSDTPS**2/ETA2-ETA*ETDDT*PSIM2-PSDTPS**2
      D23=ETDT*ETDT
      D1DDT=D1DT*(D14+D15)+D1*(D18-2.*D19+TOTHRD*D17+2.*(ALPHA2*D25
     '         /BETAO2+EO*EDDOT)/BETAO2)
      XNTRDT=XNDT*(2.*TOTHRD*D17+3.*
     1  (D25+EO*EDDOT)/ALPHA2-6.*ALDTAL**2 +
     1  4.*D18-7.*D19 )   +
     1  C1DTC1*XNDDT+C1*(C1DTC1*D16+
     1  D9*ETDDT+D10*EDDOT+D23*(6.+30.*EETA+68.*EOSQ)+
     1  ETDT*EDOT*(40.+30.*
     '  ETA2+272.*EETA)+D25*(17.+68.*ETA2) +
     1    B1*(D1DDT*D2+2.*D1DT*D2DT+D1*(ETDDT*D11+D23*(72.+54.*ETA2))) +
     1    B2*(D1DDT*D3+2.*D1DT*D3DT+D1*(ETDDT*D12+D23*(30.+30.*ETA2))) *
     1    COS2G+
     1      B3*((D5DT*D14+D5*(D18-2.*D19)) *
     1 D4+2.*D4DT*D5DT+D5*(ETDDT*D13+22.5*ETA*D23)) *SING+XGDT1*
     1         ((7.*D20+4.*EO*EDOT/BETAO2)*
     '         (C5*COSG-2.*C4*SIN2G)
     '         +((2.*C5DT*COSG-4.*C4DT*SIN2G)-XGDT1*(C5*SING+4.*
     '         C4*COS2G))))
      TMNDDT=XNDDT*1.E9
      TEMP=TMNDDT**2-XNDT*1.E18*XNTRDT
      PP=(TEMP+TMNDDT**2)/TEMP
      GAMMA=-XNTRDT/(XNDDT*(PP-2.))
      XND=XNDT/(PP*GAMMA)
      QQ=1.-EDDOT/(EDOT*GAMMA)
      ED=EDOT/(QQ*GAMMA)
      OVGPP=1./(GAMMA*(PP+1.))
      GO TO 70
   50 ISIMP=1
      EDOT=-TOTHRD*XNDTN*(1.-EO)
   70 IFLAG=0

*      UPDATE FOR SECULAR GRAVITY AND ATMOSPHERIC DRAG

  100 XMAM=FMOD2P(XMO+XLLDOT*TSINCE)
      OMGASM=OMEGAO+OMGDT*TSINCE
      XNODES=XNODEO+XNODOT*TSINCE
      IF(ISIMP .EQ. 1) GO TO 105
      TEMP=1.-GAMMA*TSINCE
      TEMP1=TEMP**PP
      XN=XNODP+XND*(1.-TEMP1)
      EM=EO+ED*(1.-TEMP**QQ)
      Z1=XND*(TSINCE+OVGPP*(TEMP*TEMP1-1.))
      GO TO 108
  105 XN=XNODP+XNDT*TSINCE
      EM=EO+EDOT*TSINCE
      Z1=.5*XNDT*TSINCE*TSINCE
  108 Z7=3.5*TOTHRD*Z1/XNODP
      XMAM=FMOD2P(XMAM+Z1+Z7*XMDT1)
      OMGASM=OMGASM+Z7*XGDT1
      XNODES=XNODES+Z7*XHDT1

*      SOLVE KEPLERS EQUATION

      ZC2=XMAM+EM*SIN(XMAM)*(1.+EM*COS(XMAM))
      DO 130 I=1,10
      SINE=SIN(ZC2)
      COSE=COS(ZC2)
      ZC5=1./(1.-EM*COSE)
      CAPE=(XMAM+EM*SINE-ZC2)*
     1   ZC5+ZC2
      IF(ABS(CAPE-ZC2) .LE. E6A) GO TO 140
  130 ZC2=CAPE

*      SHORT PERIOD PRELIMINARY QUANTITIES

  140 AM=(XKE/XN)**TOTHRD
      BETA2M=1.-EM*EM
      SINOS=SIN(OMGASM)
      COSOS=COS(OMGASM)
      AXNM=EM*COSOS
      AYNM=EM*SINOS
      PM=AM*BETA2M
      G1=1./PM
      G2=.5*CK2*G1
      G3=G2*G1
      BETA=SQRT(BETA2M)
      G4=.25*A3COF*SINI
      G5=.25*A3COF*G1
      SNF=BETA*SINE*ZC5
      CSF=(COSE-EM)*ZC5
      FM=ACTAN(SNF,CSF)
      SNFG=SNF*COSOS+CSF*SINOS
      CSFG=CSF*COSOS-SNF*SINOS
      SN2F2G=2.*SNFG*CSFG
      CS2F2G=2.*CSFG**2-1.
      ECOSF=EM*CSF
      G10=FM-XMAM+EM*SNF
      RM=PM/(1.+ECOSF)
      AOVR=AM/RM
      G13=XN*AOVR
      G14=-G13*AOVR
      DR=G2*(UNMTH2*CS2F2G-3.*TTHMUN)-G4*SNFG
      DIWC=3.*G3*SINI*CS2F2G-G5*AYNM
      DI=DIWC*COSI

*      UPDATE FOR SHORT PERIOD PERIODICS

      SNI2DU=SINIO2*(
     1   G3*(.5*(1.-7.*THETA2)*SN2F2G-3.*UNM5TH*G10)-G5*SINI*CSFG*(2.+
     2         ECOSF))-.5*G5*THETA2*AXNM/COSIO2
      XLAMB=FM+OMGASM+XNODES+G3*(.5*(1.+6.*COSI-7.*THETA2)*SN2F2G-3.*
     1      (UNM5TH+2.*COSI)*G10)+G5*SINI*(COSI*AXNM/(1.+COSI)-(2.
     2      +ECOSF)*CSFG)
      Y4=SINIO2*SNFG+CSFG*SNI2DU+.5*SNFG*COSIO2*DI
      Y5=SINIO2*CSFG-SNFG*SNI2DU+.5*CSFG*COSIO2*DI
      R=RM+DR
      RDOT=XN*AM*EM*SNF/BETA+G14*(2.*G2*UNMTH2*SN2F2G+G4*CSFG)
      RVDOT=XN*AM**2*BETA/RM+
     1      G14*DR+AM*G13*SINI*DIWC

*      ORIENTATION VECTORS

      SNLAMB=SIN(XLAMB)
      CSLAMB=COS(XLAMB)
      TEMP=2.*(Y5*SNLAMB-Y4*CSLAMB)
      UX=Y4*TEMP+CSLAMB
      VX=Y5*TEMP-SNLAMB
      TEMP=2.*(Y5*CSLAMB+Y4*SNLAMB)
      UY=-Y4*TEMP+SNLAMB
      VY=-Y5*TEMP+CSLAMB
      TEMP=2.*SQRT(1.-Y4*Y4-Y5*Y5)
      UZ=Y4*TEMP
      VZ=Y5*TEMP

*      POSITION AND VELOCITY

      X=R*UX
      Y=R*UY
      Z=R*UZ
      XDOT=RDOT*UX+RVDOT*VX
      YDOT=RDOT*UY+RVDOT*VY
      ZDOT=RDOT*UZ+RVDOT*VZ

      RETURN
      END
\end{verbatim}
-- SDP8.FOR --
\begin{verbatim}
*      SDP8                                              14 NOV 80
      SUBROUTINE SDP8(IFLAG,TSINCE)
      COMMON/E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,
     1           XNDD6O,BSTAR,X,Y,Z,XDOT,YDOT,ZDOT,EPOCH,DS50
      COMMON/C1/CK2,CK4,E6A,QOMS2T,S,TOTHRD,
     1           XJ3,XKE,XKMPER,XMNPDA,AE
      DOUBLE PRECISION EPOCH, DS50
      DATA RHO/.15696615/

      IF (IFLAG .EQ. 0) GO TO 100

*      RECOVER ORIGINAL MEAN MOTION (XNODP) AND SEMIMAJOR AXIS (AODP)
*      FROM INPUT ELEMENTS --------- CALCULATE BALLISTIC COEFFICIENT
*      (B TERM) FROM INPUT B* DRAG TERM

      A1=(XKE/XNO)**TOTHRD
      COSI=COS(XINCL)
      THETA2=COSI*COSI
      TTHMUN=3.*THETA2-1.
      EOSQ=EO*EO
      BETAO2=1.-EOSQ
      BETAO=SQRT(BETAO2)
      DEL1=1.5*CK2*TTHMUN/(A1*A1*BETAO*BETAO2)
      AO=A1*(1.-DEL1*(.5*TOTHRD+DEL1*(1.+134./81.*DEL1)))
      DELO=1.5*CK2*TTHMUN/(AO*AO*BETAO*BETAO2)
      AODP=AO/(1.-DELO)
      XNODP=XNO/(1.+DELO)
      B=2.*BSTAR/RHO

*      INITIALIZATION

      PO=AODP*BETAO2
      POM2=1./(PO*PO)
      SINI=SIN(XINCL)
      SING=SIN(OMEGAO)
      COSG=COS(OMEGAO)
      TEMP=.5*XINCL
      SINIO2=SIN(TEMP)
      COSIO2=COS(TEMP)
      THETA4=THETA2**2
      UNM5TH=1.-5.*THETA2
      UNMTH2=1.-THETA2
      A3COF=-XJ3/CK2*AE**3
      PARDT1=3.*CK2*POM2*XNODP
      PARDT2=PARDT1*CK2*POM2
      PARDT4=1.25*CK4*POM2*POM2*XNODP
      XMDT1=.5*PARDT1*BETAO*TTHMUN
      XGDT1=-.5*PARDT1*UNM5TH
      XHDT1=-PARDT1*COSI
      XLLDOT=XNODP+XMDT1+
     2           .0625*PARDT2*BETAO*(13.-78.*THETA2+137.*THETA4)
      OMGDT=XGDT1+
     1      .0625*PARDT2*(7.-114.*THETA2+395.*THETA4)+PARDT4*(3.-36.*
     2         THETA2+49.*THETA4)
      XNODOT=XHDT1+
     1       (.5*PARDT2*(4.-19.*THETA2)+2.*PARDT4*(3.-7.*THETA2))*COSI
      TSI=1./(PO-S)
      ETA=EO*S*TSI
      ETA2=ETA**2
      PSIM2=ABS(1./(1.-ETA2))
      ALPHA2=1.+EOSQ
      EETA=EO*ETA
      COS2G=2.*COSG**2-1.
      D5=TSI*PSIM2
      D1=D5/PO
      D2=12.+ETA2*(36.+4.5*ETA2)
      D3=ETA2*(15.+2.5*ETA2)
      D4=ETA*(5.+3.75*ETA2)
      B1=CK2*TTHMUN
      B2=-CK2*UNMTH2
      B3=A3COF*SINI
      C0=.5*B*RHO*QOMS2T*XNODP*AODP*TSI**4*PSIM2**3.5/SQRT(ALPHA2)
      C1=1.5*XNODP*ALPHA2**2*C0
      C4=D1*D3*B2
      C5=D5*D4*B3
      XNDT=C1*(
     1  (2.+ETA2*(3.+34.*EOSQ)+5.*EETA*(4.+ETA2)+8.5*EOSQ)+
     1   D1*D2*B1+   C4*COS2G+C5*SING)
      XNDTN=XNDT/XNODP
      EDOT=-TOTHRD*XNDTN*(1.-EO)
      IFLAG=0
      CALL DPINIT(EOSQ,SINI,COSI,BETAO,AODP,THETA2,SING,COSG,
     1          BETAO2,XLLDOT,OMGDT,XNODOT,XNODP)

*      UPDATE FOR SECULAR GRAVITY AND ATMOSPHERIC DRAG

  100 Z1=.5*XNDT*TSINCE*TSINCE
      Z7=3.5*TOTHRD*Z1/XNODP
      XMAMDF=XMO+XLLDOT*TSINCE
      OMGASM=OMEGAO+OMGDT*TSINCE+Z7*XGDT1
      XNODES=XNODEO+XNODOT*TSINCE+Z7*XHDT1
      XN=XNODP
      CALL DPSEC(XMAMDF,OMGASM,XNODES,EM,XINC,XN,TSINCE)
      XN=XN+XNDT*TSINCE
      EM=EM+EDOT*TSINCE
      XMAM=XMAMDF+Z1+Z7*XMDT1
      CALL DPPER(EM,XINC,OMGASM,XNODES,XMAM)
      XMAM=FMOD2P(XMAM)

*      SOLVE KEPLERS EQUATION

      ZC2=XMAM+EM*SIN(XMAM)*(1.+EM*COS(XMAM))
      DO 130 I=1,10
      SINE=SIN(ZC2)
      COSE=COS(ZC2)
      ZC5=1./(1.-EM*COSE)
      CAPE=(XMAM+EM*SINE-ZC2)*
     1   ZC5+ZC2
      IF(ABS(CAPE-ZC2) .LE. E6A) GO TO 140
  130 ZC2=CAPE

*      SHORT PERIOD PRELIMINARY QUANTITIES

  140 AM=(XKE/XN)**TOTHRD
      BETA2M=1.-EM*EM
      SINOS=SIN(OMGASM)
      COSOS=COS(OMGASM)
      AXNM=EM*COSOS
      AYNM=EM*SINOS
      PM=AM*BETA2M
      G1=1./PM
      G2=.5*CK2*G1
      G3=G2*G1
      BETA=SQRT(BETA2M)
      G4=.25*A3COF*SINI
      G5=.25*A3COF*G1
      SNF=BETA*SINE*ZC5
      CSF=(COSE-EM)*ZC5
      FM=ACTAN(SNF,CSF)
      SNFG=SNF*COSOS+CSF*SINOS
      CSFG=CSF*COSOS-SNF*SINOS
      SN2F2G=2.*SNFG*CSFG
      CS2F2G=2.*CSFG**2-1.
      ECOSF=EM*CSF
      G10=FM-XMAM+EM*SNF
      RM=PM/(1.+ECOSF)
      AOVR=AM/RM
      G13=XN*AOVR
      G14=-G13*AOVR
      DR=G2*(UNMTH2*CS2F2G-3.*TTHMUN)-G4*SNFG
      DIWC=3.*G3*SINI*CS2F2G-G5*AYNM
      DI=DIWC*COSI
      SINI2=SIN(.5*XINC)

*      UPDATE FOR SHORT PERIOD PERIODICS

      SNI2DU=SINIO2*(
     1   G3*(.5*(1.-7.*THETA2)*SN2F2G-3.*UNM5TH*G10)-G5*SINI*CSFG*(2.+
     2         ECOSF))-.5*G5*THETA2*AXNM/COSIO2
      XLAMB=FM+OMGASM+XNODES+G3*(.5*(1.+6.*COSI-7.*THETA2)*SN2F2G-3.*
     1      (UNM5TH+2.*COSI)*G10)+G5*SINI*(COSI*AXNM/(1.+COSI)-(2.
     2      +ECOSF)*CSFG)
      Y4=SINI2*SNFG+CSFG*SNI2DU+.5*SNFG*COSIO2*DI
      Y5=SINI2*CSFG-SNFG*SNI2DU+.5*CSFG*COSIO2*DI
      R=RM+DR
      RDOT=XN*AM*EM*SNF/BETA+G14*(2.*G2*UNMTH2*SN2F2G+G4*CSFG)
      RVDOT=XN*AM**2*BETA/RM+
     1      G14*DR+AM*G13*SINI*DIWC

*      ORIENTATION VECTORS

      SNLAMB=SIN(XLAMB)
      CSLAMB=COS(XLAMB)
      TEMP=2.*(Y5*SNLAMB-Y4*CSLAMB)
      UX=Y4*TEMP+CSLAMB
      VX=Y5*TEMP-SNLAMB
      TEMP=2.*(Y5*CSLAMB+Y4*SNLAMB)
      UY=-Y4*TEMP+SNLAMB
      VY=-Y5*TEMP+CSLAMB
      TEMP=2.*SQRT(1.-Y4*Y4-Y5*Y5)
      UZ=Y4*TEMP
      VZ=Y5*TEMP

*      POSITION AND VELOCITY

      X=R*UX
      Y=R*UY
      Z=R*UZ
      XDOT=RDOT*UX+RVDOT*VX
      YDOT=RDOT*UY+RVDOT*VY
      ZDOT=RDOT*UZ+RVDOT*VZ

      RETURN
      END
\end{verbatim}
-- THETAG.FOR --
\begin{verbatim}
      FUNCTION THETAG(EP)
      COMMON /E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,XNDD6O,BSTAR,
     1 X,Y,Z,XDOT,YDOT,ZDOT,EPOCH,DS50
      DOUBLE PRECISION EPOCH,D,THETA,TWOPI,YR,TEMP,EP,DS50
      TWOPI=6.28318530717959D0
      YR=(EP+2.D-7)*1.D-3
      JY=YR
      YR=JY
      D=EP-YR*1.D3
      IF(JY.LT.10) JY=JY+80
      N=(JY-69)/4
      IF(JY.LT.70) N=(JY-72)/4
      DS50=7305.D0 + 365.D0*(JY-70) +N + D
      THETA=1.72944494D0 + 6.3003880987D0*DS50
      TEMP=THETA/TWOPI
      I=TEMP
      TEMP=I
      THETAG=THETA-TEMP*TWOPI
      IF(THETAG.LT.0.D0) THETAG=THETAG+TWOPI
      RETURN
      END
\end{verbatim}
